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Abstract. Based on a multi-agent model, we investigate how target waves emerge 

from a population dynamics with cyclical interactions among three species. We show 

that the periodically injecting source in a small central area can generate target waves 

|>y in a two-dimensional lattice system. By detecting the temporal period of species' 

i-G \ concentration at the central area, three modes of target waves can be distinguished. 

Those different modes result from the competition between local and global oscillations 
induced by cyclical interactions: Mode A corresponds to a synchronization of local and 
global oscillations, Mode B results from an intermittent synchronization, and Mode C 
f - ^ ■ corresponds to the case when the frequency of the local oscillation is much higher than 

&\ , that of the global oscillation. This work provides insights into pattern formation in 

Cy ' biologic and ecologic systems that are totally different from the extensively studied 

. ■ diffusion systems driven by chemical reactions. 
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1. Introduction 

Spatially distributed excitable systems are widely investigated owning to their biological 
significance of long-rang signal transmission through self-sustained waves [H [21 El El El E] . 
Pattern formation of wave propagation in excitable systems have been well studied 
[7j. Especially, target waves induced by noise have been found in the recurrent single 
species' population dynamics [6]. However, the formation mechanism of target waves 
in populations dynamics driven by multi-species' competitive interactions are not very 
clear. 

Competitive interactions in natural and social systems consisted of many elements 
(e.g., various biological species, political parties, businesses, coupled reaction chemical 
components, bacterial production bacteria, etc.) play an important role in evolutionary 
processes and pattern formations P, El [10, [HI H21 EB], and lead to the emergence of 
spatial patterns. In particular, target waves are commonly observed in those systems 
[HI EEEl ESI El EEE]- Previous works reveal that patterns could emerge in reaction- 
diffusion systems if one of the substances diffuses much faster than others [191 [201 [21]. F° r 
example, the bromous acid diffuses much faster than ferrion in the Belousev-Zhabotinsky 
(BZ) reaction and the cyclic AMP diffuses much faster than membrane receptor in the 
dictyostelium discoideum. However, many pattern formations of mobile population in 
ecosystems, such as migrating animals and running bacteria, can not be explained by the 
react ion- diffusion mechanism, since the diffusion speed induced by individual mobility is 
the same for all substances [221 [23]. Furthermore, partial differential equations (PDEs) 
have been proposed to describe the evolution of pattern formation in reaction-diffusion 
systems, such as the Oregonator model for BZ reaction [21]. Based on those PDEs, one 
can analyze the stability of patterns by using the mean field theory, but the reliable 
detailed information is quite limited. Therefore, it is necessary to use an agent-based 
models to describe the pattern formation of mobile population [H] . 

The cyclic predator-prey model provides a terse description of competition among 
population of different species. Experimental study [22] has revealed that the 
mechanism of rack-paper-scissors game can promote the biodiversity of three strains 
of Escherichia coli. Very recently, Reichenbach et al proposed a rack-paper-scissors 
game of mobile populations [8], where the individual mobility displays a critical effect 
on species diversity. When mobility is below a certain value, all species coexist and 
form spiral waves. In contrast, above this threshold biodiversity is jeopardized. Indeed, 
previous investigations show that the cyclical competition mechanism and low mobility 
of individuals maintain the biodivesity in ecosystems [HI [25], [22], [26] [27] . In this paper, 
we investigate the pattern formation based on a cyclic predator-prey model with mobile 
individuals. Target waves are observed from the recurrent dynamics driven by three 
species' competitive interactions. As a result of global oscillation, a transition from 
disordered state to an ordered spacial structure occurs with the increasing of injection 
period in the vortex of target waves. By detecting the temporal period of species' 
concentration at the central area, three modes of target waves can be distinguished. 
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Those different modes result from the competition between local and global oscillations 
induced by cyclical interactions: Mode A corresponds to a synchronization of local 
and global oscillations, Mode B results from an intermittent synchronization, and 
Mode C corresponds to the case when the frequency of the local oscillation is much 
higher than that of the global oscillation. Though the present model only concerns 
the microscopic interactions among individuals, target waves emerge in a macroscopic 
level. Furthermore, our work provides a feasible method for pattern control and pattern 
selection. 

2. Stochastic cyclic predator-prey model 

Based on the previous works of Reichenbach et al [El|25], we introduce a cyclic predator- 
prey model as follow: Mobile individuals of three species (marked by 1, 2, 3) locate in 
the nodes of a two-dimensional regular lattice with no flux boundary. Here a node stands 
for an 1 x 1 square. Each node can be occupied by at most one individual of a species or 
a vacancy (denoted by V) representing resource. There are three kinds of interactions, 
namely exchange, predation, and reproduction, which only occur between neighboring 
nodes. Exchange. — An individual could exchange positions with one of its neighbors 
at a rate a due to its mobility. Predation. — 1 beats 2 and 2 becomes a vacancy at a 
predation intensity /3, in the same way, 2 beats 3, and 3 beats 1. Reproduction. — An 
individual can reproduce an offspring to a neighboring V node at a rate 7. 

Unlike the deterministic approach which regards the time evolution as a continuous 
process, here, the applied stochastic approach regards the time evolution as a kind of 
random-walk process. In this model, reactions occur in a random manner: exchange 
happens with probability a /(a + (3 + 7), predation with probability ft /(a + j3 + 7), and 
reproduction with probability j/(a+0 + j). We use a standard algorithm for stochastic 
simulation, which was developed by Gillespie [28], [29] . 

3. Results 

Inspired by the experiments of growth of Escherichia coli on a petri dish [15] , we apply 
the dynamical model of three populations' cyclical interactions [8]. In an L x L lattice, 
each node presents a space occupied by a vacancy or an individual of a species, and 
at each time step a randomly chosen individual interacts with one of its four nearest 
neighbors, which is also determined randomly according to the corresponding predation 
intensity. The evolving process is implemented by using the Gillespie algorithm [281 ES] 
in which one Monte Carlo (MC) time is defined as the time period during which all 
the individuals have been chosen once on average. Hereinafter, the resolution of time is 
defined as one MC time. Initially, the lattice is wholly occupied by vacancies, and then 
populations of a species are injected periodically in a central area, namely the injected 
area. An illustration of the injected area is shown in Fig. [U In this paper, the injected 
radius, R, is set as 10.5 and the lattice size, L, is 1024. The injected species revolve in 
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Figure 1. Illustration of the injected area in an 11 x 11 lattice. Here, R is set as 3.5, 
and the lattice completely inside the circle, marked by the blue color, is the injected 
area. In the simulation of this paper, the injected radius R is fixed as 10.5. 



the order 1, 2, 3, 1, 2, 3, ... For example, populations of 1 are injected at t = 0, then 
2 at t = T , 3 at t — 2T , and 1 at t = 3T , etc. In this case, the period of injecting is 
T in = 3T , where T is the interval time between injections of two species. 

Mobility is an important character for most population dynamics such as bacteria 
swimming and tumbling. The mobility in the present model corresponds to the action 
that an individual moves to a neighboring empty node (V node), which is defined as: 



M = 2a/N, 



(1) 



where N = L 2 is the number of nodes in the system. In this paper, L = 1024, (3 = 1 
and 7 = 1 are fixed. 

With no-flux boundary conditions, extensive computer simulations are performed 
and the system displays abundance phenomena of pattern formation with various 
injection periods. In order to describe the spatial structure of pattern formation, 
we define p = Nq/N, where Nq is the number of individuals in the largest adjacent 
population. Here, an adjacent population is a set of individuals of one species where any 
two individuals can be connected through a sequence of nearest neighboring relations 
within this set. Apparently, patterns break when p approaches to 0, while one species 
occupies the whole system when p reaches 1. As shown in Fig. [2]^d), given M = 10 -5 , 
with the increasing of To, the system undergoes a transition from a state of coexistence of 
three species to a state occupied by one species. The pattern formations corresponding 
to the points a, b, and c in Fig. EJ^d) are shown in Fig. Eta), E](b), and[2](c) respectively. 
When T is small (T < 300) the three species form many small spiral waves which can be 
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Figure 2. The panels (a), (b), (c), (e), (f) and (g) display different pattern formations 
generated in the present model. The panels (d) and (h) report p as a function of T . 
(a), (b), (c), (d) for M = 1CT 5 and (e), (f), (g), (h) for M = 1CT 4 . (a) and (c) for 
T = 300, (b) and (f) for T = 500, (c) and (g) for T = 1000. The cases presented 
in (a), (b), and (c) are marked by by a, &, and c in (d), and the cases presented in 
(e), (f), and (g) are marked by e, /, and g in (/i). In (d), the relation p(Tq) can be 
approximately linearly fitted in the interval To £ (500, 1800). In (g), y is denotes the 
spatial distance between neighboring wave fronts for the same species. 



considered as a disordered state in global level. Fig. [2(b) and Fig. [2(c) show target waves 
emerging for To = 500 and To = 1000. The system becomes more and more ordered, 
reaches a consensus state when T > 2000. In contrast, for M = 10 -4 , Fig. [2(h) shows an 
unexpected decline with the increasing of To, and the pattern formations corresponding 
to the points e, /, g, are shown in Fig. 12(e), 12(f), and 12(g), respectively. 

It is worth mentioning that all patterns in our results are robust after the systems 
have reached steady states. It seems that target waves emerge because of the increasing 
of T in . A question is naturally taken into consideration: How T in affects the pattern 
formation. Furthermore, one may ask why p is monotonously increasing with the 
increasing of T in Fig. [2(d), while non-monotonously varying in Fig. [2(h). We calculate 
the proportion of species 1 in the injected area. With T = 600, this proportion does not 
display periodicity for M = 10~ 6 (see Fig. [3(a)), while varies periodically for M = 10 -5 
(see Fig. [3(b)) as well as for M = 10~ 3 (see Fig. [3(d)), and shows quasi-periodical 
behavior for M = 10~ 4 (see Fig. [3(c)). The period is T = 2T in = 3600 for M = 10" 5 , 
and T = T in = 1800 for M = 10" 3 . Fig. [3(e), [3(f), [3(g), and[3(h) show (the varying of) 
three species' proportions in the whole system one-to-one corresponding to the evolving 
processes of Fig. [3(a), [3(b), [3(c), and [3(d). Clearly, Fig. E(a)-[3(d) show the local 
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Figure 3. The panels (a)-(d) report the proportion of species 1 in the injected area 
for different mobilities, while (e)-(h) display proportions of three species in the whole 
system. Panels (e)-(h) are one-to-one corresponding to panels (a)-(d). To is fixed 
as 600, M = 10~ 6 for (a) and (c), M = 10~ 5 for (b) and (f), M = 10~ 4 for (c) 
and (g), M = 10~ 3 for (d) and (h). The panel (k) shows different system states from 
A/ = 10~ 10 to M = 10 -2 , given To = 600. The insets (i) and (j) display the comparison 
of species l's proportions in the injected area (black) and the whole system (red),(i) 
for M = 10~ 5 (corresponding to (b) and (f)) and (j) for M = 10 -3 (corresponding to 
(d) and (h)). 



oscillations in the injected area, while Fig. [3]^e)-mh) display the global oscillations. 
Fig. [3]^i) and Fig. [3](j), corresponding to the cases shown in[3]^b) and[3^d), compared the 
local and global oscillations of species 1. In Fig. E](i), the local and global oscillations 
have the same period, while in Fig. H2(j), the frequency of the local oscillation is much 
higher than that of the global oscillation. We say, in the case shown in Fig. [3](b), [3^f) 
and[3](i), the local and global oscillations are synchronous at T = 2Tj n , while in the case 
shown in Fig. [3]^d), [3](h) andEfj), they are non-synchronous. In between, the system 
displays quasi-periodical behavior, and the local and global oscillations are intermittent 
synchronous (see, for example, Fig. [3tc) and^g)). According to the above discussion, a 
qualitative classification of system states versus mobility is shown in Fig.[3](k), including 
non-period, period T = 2T in , quasi-period, and period T = T in . 
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Figure 4. The panels (a)-(d) show pattern formations corresponding to the cases 
shown in Fig. EJa), (b), (c), and (d) respectively. 



As shown in Fig. HI there are four qualitatively different states, and the patterns in 
Fig. 1Kb), Fig. BJc) and Fig. |D[d) are named as Mode A, Mode B and Mode C. When 
a species is injected, individuals of this species will interact with other species around, 
(i) For extremely small value of M (< 3.0 x 10 -6 ), individuals out of the injected 
area distribute almost homogeneously due to a relatively high rate of predation and 
reproduction. The injected individuals, surrounded by other two species, are segmented 
into fragments. Therefore, as shown in Fig. Ufa) the local oscillations have no period, 
and no global ordered pattern formats, (ii) For small value of M (< 10~ 4 ), individuals 
out of the injected area are inhomogeneously distributed due to a relatively low rate 
of predation and reproduction. The injected individuals are surrounded by only the 
prey species, and each species takes up the injected area for 2T in /3, thus the period of 
local oscillations is 2T in . The boundary between the two species near the injected area 
moves from the center to the whole system, inducing a global oscillation with period 
2T in too. This case is shown in Fig. H(b) and referred to as Mode A. (Hi) For medium 
value of M (10 -4 < M < 5.0 x 10~ 4 ), although individuals out of the injected area are 
inhomogeneously distributed due to a relatively low rate of predation and reproduction, 
the injected individuals are surrounded by all three species. No species takes up the 
whole injected area and the period of local oscillations is 2T in sometimes (see Fig. El(c)). 
In other words, local oscillations are confined in the injected area sometimes. Therefore, 
the local oscillations is quasi-periodic and intermittently synchronized with the global 
oscillations. This case is shown in Fig. H](c) and referred to as Mode B. (iv) When M is 
sufficiently large (> 5.0 x 10~ 4 ), the relative rate of predation and reproduction is very 
low. Three species coexist near the boundary of the injected area. Because the number 
of injected individuals is quickly depressed, no species dominates the whole injected area 
(see Fig. EJ^d)). Local oscillations are confined in the injected area, and individuals of a 
species will aggregate near the boundary of the injected area, and then move from the 
center to the whole system. Since the aggregation and propagation takes a long time, 
the period of global oscillations is much longer than that of local oscillations. Target 
waves formatting in this manner are called Mode C. 

The wavelength is defined as A = y/L, where y is denotes the spatial distance 
between neighboring wave fronts with the same species in the lattice as shown in 
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Figure 5. (a) Wavelength of target waves as a function of M for different values of 
injection period Tq. The inset shows A ~ logM 1 / 2 . (b) Different system states in 
A I - T Q space. 



Fig. |2(g). Comparing Fig. H(b), Fig. IH(c), and Fig. H(d), one can find that wavelengths 
of target waves in both Mode A and Mode C are longer than that in Mode B, and 
wavelengths of target waves in Mode B are inequable. It can explain why Fig. E](h) 
have a trop. The reason is target waves belong to Mode B for T > 300 and T < 800 
at M = 10 -4 , while others belong to Mode A or Mode C, and wavelengths of Mode B 
are smaller than that of Mode A and Mode C. Moreover, the target waves become 
unstable when T goes across the critical point distinguishing two different modes, and 
may evolve into spiral waves which is also observed in a system of a single species, 
Dictyostelium system [6| [30]. 

As shown in Fig. 0(a), wavelengths of target waves display non-monotonic behavior 
with the mobility of individuals for T = 600 and T = 400, while increase with M 
in terms of A ~ log M 1 ' 2 for To=300. This phenomena are different from the cases of 
spiral waves [8], in which A ~ M 1//2 . The non-monotonic increasing with the mobility 
T = 400 and T = 600 is caused by the target waves of Mode B (see red blocks in 
Fig.EKb)). 



4. Summary 

Cyclical competition, an important principle of ecology and society, leads to complex 
phenomena of maintenance of ecological bio-diversity and emergence of cooperation in 
society, where spatial distribution of populations plays a critical role [HH [32} [331 134] . 
It has been found that multi-species' competition induces complex spatiotemporal 
structures in ecological systems. Especially, the 3-species evolving game, namely rock- 
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paper-scissors game or cyclic predator-prey model, reveals spiral wave structures 
However, target waves are rarely reported. In the current model, target waves emerge 
under a suitable initial condition, which is in accordance with the observations in 
real bacteria systems [61 [30]. Although target waves were also observed in many 
excitable media, in our model, individuals of three species move at the same rate, 
which distinguishes our model from the diffusion systems driven by chemical reactions 
where different substances have different moving speeds. Indeed, our results suggest 
that target waves can be driven by the competition of local and global oscillations. This 
newly reported mechanism could provide insights into target waves formation in biologic 
and ecologic systems. 

In addition, our results can be applied in pattern control. Actually, the system 
makes stochastic resonance to the periodical injection fluid in Mode A, which can be 
used to control the global behaviors since the injection period is controllable. Recently, 
the periodic injection method is also used in complex Ginzburg-Landau system for 
the purpose of spatiotemporal chaos control [33 EH]- The periodical injection for 
pattern control has great significance in potential applications of toothful medication 
transportation. 
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